//---   The programme "NonSinusoidalDrivenDampedOscillator", Ebrahim Foulaadvand, 08 August 2012 -----

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;


double a (double &t,double &x,double &v); // a is the accelration function

double  tau=0.01,T=100,ac,k=9.,m=1.,w0=sqrt(k/m),x0=0.,v0=0.,gamma=0.5,beta=gamma/(2*m),omega=2.,
        period=2*M_PI/omega,delta,A0=1.,E,kx1,kx2,lx1,lx2,tmid,xmid,vmid,Anumeric=0.;

// x0=initial position ; v0=initial velocity ; tau=timestep; T=simulation time,
// gamma=damping coefficient; omega= frequency of external force


int N=(T/tau),

main()

{


vector<double> v(N+1,0),x(N+1,0);

ofstream xER("xER nonsinusoidal omega=2.0.plt");
ofstream xRK2("xRK2 nonsinusoidal.plt");
ofstream aER("aER nonsinusoidal.plt");


//------------------------ Euler-Richardson Method------------------------

x[0]=x0; v[0]=v0; Anumeric=0;

for(int t=0;t<N;t++){

tmid=t*tau + 0.5*tau;
xmid=x[t] + 0.5*v[t]*tau;
vmid=v[t] + 0.5*tau*a(t*tau,x[t],v[t]);

v[t+1]=v[t] + tau*a(tmid,xmid,vmid);
x[t+1]=x[t] + tau*vmid;


xER<<t*tau<<" "<<x[t]<<endl;
aER<<t*tau<<" "<<a(t*tau,x[t],v[t])<<endl;

if (t*tau>(10/beta) && t*tau<(10/beta) + 2*2*(M_PI/omega) && x[t]>0 && x[t+1]>x[t])  Anumeric=x[t+1];


}//end t


//cout<<"Anumeric ER= "<<Anumeric<<endl;


//------------------------- Runge-Kutta 2nd order Method ------------------------

x[0]=x0; v[0]=v0;

for(int t=0;t<N;t++){

kx1=tau*v[t];
lx1=tau*a(t*tau,x[t],v[t]);

kx2=tau*(v[t]+lx1);
lx2=tau*a(t*tau+tau,x[t]+kx1,v[t]+lx1);

x[t+1]=x[t]+0.5*(kx1+kx2);
v[t+1]=v[t]+0.5*(lx1+lx2);

xRK2<<t*tau<<" "<<x[t]<<endl;


if (t*tau>(10/beta) && t*tau<(10/beta) + 2*(M_PI/omega) && x[t]>0 && x[t+1]>x[t])  Anumeric=x[t+1];

}//end t


//cout<<"Anumeric RK2= "<<Anumeric<<endl;getch();


cout<<"end."<<endl;getch();

return 0;

}

//----------------------------------------------------------------


double a (double &t,double &x,double &v) {

if ( fmod(t,period)<=0.25*period || fmod(t,period)>=0.75*period )

{ return ac=-k*x/m-(gamma/m)*v + A0*cos(omega*t); }

if ( fmod(t,period)>0.25*period && fmod(t,period)<0.75*period )

{ return ac=-k*x/m-(gamma/m)*v; }

}


/*
double a (double &t,double &x,double &v) {

return ac=-k*x/m-(gamma/m)*v + A0*cos(omega*t);

}
*/